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Abstract 



i-^ The existing solution for the Langevin equation of an anisotropic fluid allowed the evaluation of 

^ the position dependent perpendicular and parallel diffusion coefficients, using Molecular Dynamics 

4_> data. However, the time scale of the Langevin Dynamics and Molecular Dynamics are different and 

c/3 an anzat for the persistence probability relaxation time was needed. Here we show how the solution 

C^ for the average persistence probability obtained from the Smoluchowski-Fokker-Planck equation 

1 (SE), associated to the Langevin Dynamics, scales with the corresponding Molecular Dynamics 

^ quantity. Our SE perpendicular persistence time is evaluated in terms of simple integrals over the 

I equilibrium local density. When properly scaled by the perpendicular diffusion coefficient, it gives 
a good match with that obtained from Molecular Dynamics. 
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I. INTRODUCTION 

The understanding of the dynamics of anisotropic fluids is fundamental in the study 
of many chemical and biophysical processes, particularly those occurring under nano- 
confinement. Molecular quantities like the mean square displacement (MSD), the mean first 
passage time (MFPT), and probability functions as the survival or persistence probability, 
can be readily obtained from a regular Molecular Dynamics (MD) simulation. However, 
phenomenological quantities are bounded to the equations and approximations that defined 
them. Thus, one commonly has to invoke some sort of approximation or condition for their 
validity. For instance, the diffusion constant is defined through the stochastic Langevin 
equation as a ratio of the approximate fluctuation force to the mechanical friction force in 
the so called Smoluchowski time scale. On a previous work jjj , we presented an analytical 
solution for the Langevin equation of an anisotropic fluid, next to an attractive wall, 
which allowed the evaluation of the position dependent perpendicular diffusion coefficient. 
The mean squared displacement (MSD) and the average persistence time t{z), needed to 
compute this coefficient, were obtained from a modification of the virtual layers molecular 
dynamics method (VLMD) of Liu et a/, p]. 

A serious difficulty found was that the time scale of the Langevin Dynamics and of the 
Molecular Dynamics did not match. In fact, it was discussed by Burschka et al. [3], Harris 
[1], and Razi Naqvi et al. [5J , that the conditional probability obtained as the stationary 
solution of the Fokker-Planck equation, associated to the Langevin equation for a system 
with absorbing boundaries conditions fails to vanish at those boundaries. For the same 
space sampling layers width L, the Langevin persistence probability is then known to be 
lower than the corresponding MD persistence probability. That is, MD particles reach 
the boundary faster than Langevin particles. In another words, the Langevin dynamics 
(LD) average persistence time, r^^, for a given sampling layer width, is lower than the 
corresponding molecular dynamics time r^^^ . This is a systematic inconsistency, found even 
for virtual absorbing layers located at the bulk. It has been handled in the literature by 
shifting the boundaries away. In fact, Liu, Harder and Berne [2] overcame this difficulty, 
by running simultaneous dual numerical LD and MD simulations, changing the value of L, 
until the survival probabilities from both simulations matched. Alternatively we appealed 



to the simplest ansatz of fixing the layer width and shifting the r^^^ time instead [T] . 

In this paper we take a closer look at this problem, by solving the backwards 
Smoluchowski-Fokker-Planck equation (SE) associated to the anisotropic Langevin equa- 
tion and comparing the resulting r^^ {z) with the r^,^^ {z) measured from the VLMD method. 
In the next section we briefly review the probability concepts in the VLMD method. After 
establishing the relation between the mean first passage time and the persistence time, we 
shall obtain an expression for the SE average persistence time. Finally in the last section, 
we shall discuss the scaling between the SE and MD dynamics, carrying out calculations 
for a simple Lennard- Jones dense fluid next to an attractive wall. 



II. VIRTUAL LAYER MOLECULAR DYNAMICS METHOD 

In this work we will use the VLMD method as discussed previously[T] . It is a modifica- 
tion of the one introduced by Liu et al. (LHB) [2] and Thomas et al. [6] . Briefly, once the 
dynamics is equilibrated, the total simulation time s is divided in discrete time steps As 
and partitioned in J blocks containing n^ax time steps. The MD simulation time interval 
after n steps is then t = nAs. The z coordinate perpendicular to the attractive walls is 
also divided in discrete virtual layers of width L. We then consider the set of particles 
that stay in a given layer located at Za, i.e. z(t) G {a,b} = {za,Za + L}, during the time 
interval t, spanning between the simulation time sq and s = sq + t. The initial number of 
particles in the layer at s = Sq is A^(0) = A^(so), and the number of particles in the set, still 
in the layer after the interval t, is N{t) = N{sq + t). The maximum time interval used to 
evaluate quantities inside the layers was tmax = nmaxAs. After rimax steps the algorithm is 
re-initiated to measure the dynamics in the layer, setting sq = s again. This layer sampling 
is repeated J times. We denote the number of particles in the set that stay in a layer in the 
j^^ layer sampling or repetition as Nj{t). 

In the VLMD method, the dynamical properties of the anisotropic fluid are evaluated 
layer by layer using absorbing boundary conditions. That is, particles that exit a given 
layer are not further counted, even if they reenter the region. The dynamics in the z 



direction is governed by the behavior of P{z, t; Zg), the joint probabihty of finding a particle 
at position z^ at time t = and then at position z at time t. For an anisotropic system 
this function depends on the position where the layer is located, a= z^, and on the width L 
of the space interval. The validity and accuracy of the VLMD method rely on the proper 
choice of L[ll [2] . It was found that if L is too large the mean force within the layer changes 
more than required by the method. While, if it is too small, particles tend to escape from 
the virtual layer, affecting the proper statistics of the MD. For the fluid conditions used, L 
equal to 0.5 diameters was found to be a good compromise. 

Since the movement in the z direction is assumed to be independent of the movement in 
the X and y directions, the joint probability function is given as 

P{z, t,Zg) = P{z, t\zg)pa{z,), (1) 

where P{z,t\zg) is the bayesian conditional probability that the particle is located at z at 
time t, given that it was at z^ at time t = 0, in a layer containing z^. Paiz^) is the probability 
for a particle to be at position z^ at time t = 0, normalized in the given region. It can be 
obtained directly from the local particle density p{z). 

Pai^o) = -J^- (2) 

The normalization of the joint probability is by definition the persistence probability, 
V(t,Za), in the diffusion domain Za < z < Za + L 

PZa+L rZa + L 

V{t,Za)=j dz dZgP{z,t;Zg). (3) 

J Za J Za 

It measures the average probability that after a time interval t a particle still remains 
inside a layer located at Za, independently of the initial position. V(t; Za) is often also called 
the survival probability [2] , however, we prefer to call it the persistence probability to 
distinguish from the survival probability G{t, z,^; z^), as defined in the literature [H [8] 

G{t,z,;zJ= I P{z,t\z,)dz. (4) 

J a 

G(t,ZQ;z^) is the probability that at time t, the particle, initially located at Zg, remains 
anywhere in the region {a, b}. So, while the survival probability is a function of the initial 



position, the persistence probability is independent of it. 

In a MD simulation, the probability P(t, Za) is obtained by averaging Pj{t) = Nj(t)/Nj{0) 
over the J repetitions 

As it has been shown by Molecular, Generalized Langevin and Langevin Dynamics 
V{t,z^) can be phenomenologically approximated with a high degree of accuracy by an 
exponential decay [H El E] , 

V{t,zJ = e-''^''^\ (6) 

where t{z^) is the the relaxation time of the persistence process, also referred to as the 
average persistence time. The MD t{z^) is evaluated from Eq. ([6]), by fitting the numerical 
V{t,zJ MD data. 

All other physical properties of interest are evaluated as averages over the joint prob- 
ability. For instance, the mean square displacement (MSD) in a layer is obtained as the 
average 

((^W -^o)'>, = / " dz f^ dz,[zit)-z,fPiz,t,z,). (7) 

J Za J Za 

This expression also holds for the x and y directions. In the MD simulation the mean square 
displacement (MSD) within the layer at a is evaluated summing over all the Nj{t) particles 
in the set 



m-=m\ = j^j^^ 






(8) 



III. MEAN FIRST PASSAGE TIME AND PERSISTENCE PROBABILITY 



In this section we will develop the relation between the MFPT and the mean survival time 
of particles dwelling in the layer. Using Eqs. ([I]) and Q, we rewrite Eq. ([s]) as an average 
over the initial position and find that a particle dwells in a given layer with a persistence 



probability 



V{t;za) 



Za+L 

dz 



Za+L 



n \ ? o ) 



Pa{Z„, 



dz^G{z„t)p^{z^] 



(9) 
(10) 



In the context of this note, the average persistence time that the particles spend in motion 
within the virtual layer is then given, according to Eq. (|6| by the integration of V{t,z^) 
over the time interval. Using Eqs. (|9| and (10) we get 



T(Z„ 



dtVit,z^ 



dz,. 



dtG{z^,t) 



Pai^o] 



dZo[tMFpi^o)]Paiz,] 



(11) 



where we have defined the mean time it takes a particle to reach the boundary given that 
it started at z,, as the mean first passage time tj^jpp{z^-^) (MFPT)[7j 



^MFP\^0) 



dtG{Zf^,t). 



(12) 



Note that when the virtual boundaries of the layer are absorbing, the MFPT is actually 
a mean exit time. Thus, according to Eq. (11), t(z^) is simply the space average of tj^^^p^z^^) 



^(^a) = (^MFp(^a))a- 



(13) 



The MFPT have the nice feature that they can be evaluated in terms of the equilib- 
rium particle density. In earlier work, we have associated the MFPT with effective diffusion 
constants for confined ionic and non ionic fluids, and with local diffusion coefficients in 
anisotropic nano-confined molecular dense fluids [H EHTO] . So, in the next section we will 
obtain the expressions for the MFPT from the Smoluchowki-Fokker-Planck equation asso- 
ciated to the Langevin equation. 



IV. LANGEVIN DYNAMICS MEAN FIRST PASSAGE TIME 

The z component of the anisotropic Langevin-like equation, can be written as[I], |2] 

dv{z] 



m- 



dt 



-l^Az) + F{z) + a{z)i{t). 



(14) 



It describes the dynamics of a fluid particle of mass m located at position z{t) with velocity 
v(z), in the presence of and external force F{z) and a random or fluctuation force a{z)C,{t). 
7^^ is the z component of the friction tensor of the fluid and, ^(t) is the usual (5-correlated- 
zero-mean white noise, resulting from the collisions with the rest of the fluid. The coefficient 
a{z) measures the amplitude of the fluctuation force. In terms of a potential of mean force 
W{z), felt by the particles due to the interactions with the walls of the container or and 

external field, the force is given by F{z) = ^^ = k^T '^^ , where p{z) is the local 

particle density. Here, kg and T are the Boltzmann constant and temperature, respectively. 



To write Eq. (14), we used the standard assumption that the friction tensor, T, is diagonal 
with 'jzziz) 7^ ixxiz) = 1yy{z) [T! l2] , ueglectiug the off-diagonal terms j^ziz) = 1yz{z) = 0. 
Here we used 7(2;) for the perpendicular or transverse, zz diagonal element of the friction 
matrix, 7(2;) = 7^^. It has been shown that the parallel xx and yy components are also 
dependent on the z position, but to a much lesser extent [1]. Therefore, we shall only consider 
in this paper the transverse component of the anisotropic diffusion coefficient. 

In what follows, we suppose the velocities of the particles attain a canonical distribution 
much faster than positions. Thus, at the time scale of positions, velocities attain the steady 
state, -^(iv(z)/(it = 0, so the instantaneous relaxation approximation (IRA) [15j can be 

Izz 



applied to Eq. ( 14 ) to get 

dz{t) 



v^(z) + v/2DMeW, (15) 



dt 

where we defined the space dependent diffusion coefficient D{z) = D^^, to relate the fluc- 
tuation and the dissipation forces, as D{z) = |(^^)^ and, we introduced the drift velocity 
v^ (z) = -T^y . The validity of this approximation, even for ionic fluids, was successfully tested 



by Jonsson and Wennestrom|16j . Equation (15) gives the velocity of the particle under a 



stochastic potential fleld. This is also referred to as the high friction limit approximation 



and extensibly used in the literature OIH]. Thus, if the stochastic differential equation (15) 
is interpreted in Ito's sense [7] , the corresponding forward Fokker-Planck (FP) equation 
for the conditional probability P{z,t\Zg)=P{z,t\Zg,tg = 0) of flnding the particle in position 
z = z{t) at time t, given that it started to diffuse from an initial position z^ at tp = is|14j 

^^^^ = ^ [D{z)PMz,)] - I [v^(z)P(.,t|zJ] . (16) 

When one postulates that, at equilibrium, the local particle density, written as a Boltz- 
mann distribution p{z) = Pg exp{—(3W{z)), is to be a stationary solution of the FP equation, 



Eq. (16), and that there is no net flux of particles, one gets the fluctuation- dissipation the- 



orem in the form 



dD{z) 
dz 



PF{z) 



^-^ Diz) 



(17) 
(18) 

where v^(2:) = f3F(z)D(z), and as above, v^(z) = f3F{z) -^ with (3 = 1/kgT. Since 
we are interested in anisotropic systems, where —^ is nonzero, in general v^{z) ^ v^(z). 
Therefore, the Sutherland- Einstein relationship, D{z) 



-^, does not apply. Instead, 



substituting Eq. (18) in the FP equation, gives the so called Smoluchowski-Fokker-Planck 



(SEP) equation in terms of D{z) and v^{z), which can be written as a continuity equation 
in the form 



dP{z,t\z^) dj{z,t\z,) 



dt 



dz 



with the flux j{z,t\Zf^) given as 



j{z,t\Zo) =v^{z)P{z,t\z,) - D{z) 



dP{z,t\z,) 
dz 



(19) 



(20) 



This is also commonly referred to as the forward Smoluchowski equation (SE). Particle 
dynamics and fluctuations occurring in biological and liquid environments are often 
described well by this equation[l7] . The main objective of this paper is to establish the 
time scale associated to this equation, when compared to MD. 



Defining the Smoluchowski operator 



C{z) 



d 
dz 



p{z)D{z] 



d 1 
dz p{z) 



and its adjoint 



Ciz,^ 



p{zjdz, dz, 

one can rewrite the previous forward SE as 

dPiz,t\z,) 



dt 



Ciz)Piz,t\z,), 



and the corresponding adjoint or backward SE 

dP{z,t\z,) 



dt 



C*{z,,)P{z,t\z,] 



(21) 



(22) 



(23) 



(24) 
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where the operator C*{z^^) acts on the initial position z^. In the problem of first passage 
times, fully described elsewhere [3 |H] , the objective is to know how long a particle, whose 
position is described by the above equations, remains in the region {a, 6}. 

From the definition given in Eq. Q, a simple integration of Eq. (24) over 2;, gives a 
differential equation for the the survival probability 

^^|^ = £*(^JG(z„,t), (25) 

which must to be solved with the initial condition G{z^^,t = 0) = 1. Carrying out the same 



integration over both sides of the continuity equation, Eq. (19), we get 

^^^=jia,t\z,)-jib,t\z,). (26) 

The right hand side gives the flow out of the layer. From this, it can be seen that the 
first passage process to the absorbing layers is in fact an exit process. 



The mean exit time or mean first passage time is then evaluated, from Eq. (12), by inte- 



grating the evolution equation for G{zg,t), Eq. (25). So, the MFPT satisfies the differential 
equation 

C*{z,)t,,,,{z,) = -l. (27) 

The boundary conditions for these equations must fulfill the requirement that the 
virtual layers are absorbing at a and b, in the sense that particles are counted out 
once they reach the boundary. While some difficulty might exist when dealing with 
the probability function P{z,t\zQ) in the SE [2j and even in the stochastic free diffusion 
equation [H El US] , the absorbing boundary condition for a survival or exit probability 
is straightforward, namely, G{a, t) = G{b, t) = 0, and for the MFPT, t,^pp{a) = t^j^pip) = 0. 

The solution of the differential equation for the MFPT with absorbing boundary condi- 
tions follows 

, , , (r° ^) /: ^ /: p(-)^- ' (ii: jl?) £•' 4g r ^(->^- „,, 

'^MFpi^o) - Jb~Yz ■ ^ ^ 

Ja VC^) 

Here, the potential of mean force was related to the local particle density p{z^^) and we 
defined the extensive diffusion density function ipi^z^) = p{zq)D{z^). 



Following the VLMD approximations, we now assume that the layer is small enough 
that the property D{z) is constant for z G {a,b}. That is, we assume that the Sutherland- 
Einstein relationship holds locally in the virtual layer, so that -^^ = D{z) = D{z^), where 
the constant D{z^) is the diffusion coefficient evaluated at the left boundary of the virtual 
layer at position z^. 



From Eq. (13), we identify the relaxation parameter r„^ of the persistence probability, as 



the position-averaged MFPT. For an anisotropic system this position average depends on 
the position a = z^ where the region is located and, on the width L of this space interval. 
We finally get 



^SEi^a) 



where the integral H[y, s] is 



:'^)Hi.„b]- 



p(y)J 



.lsi))*l"-^.l 



£>('-.) />('-)* 



b _dz_ 
a p{z) 



(29) 



H[y,s] 



dz 



p{x)dx. 



(30) 



ly P{z) J a 

Analytical results can be drawn from the previous equations whenever the functional 
form of the local particle density p{z) is known. A limiting special case is obtained for a 
region located on the bulk reservoir away from the walls, where the potential of mean force 
vanishes. Then, using p{z^) = p^^^^, in Eq. (28) and Pa(-2()) = 1/i^ in Eqs. (11) and (29) we 
get 



&^'' (zA 

MFP ^ J 



Jjulk 



{z^ -a){h~z^) 



2f)bulk 



L' 



(31) 
(32) 



SE '^2£)bulk ' 

where D^^'-^ = /}„ is the homogeneous bulk fluid diffusion constant. It is easily evaluated, 
for example, from the long time limit of the MSD in a regular MD simulation. 

V. RESULTS AND DISCUSSION 

We studied the same anisotropic system as in reference [1], namely a dense Lennard- 
Jones fluid with a bulk reduced density of pbuik = 0.69, using 1296 molecules, at a reduced 
temperature of T* = ksT /epF = 0.75, next to a highly interacting 9-3 LJ smeared wall 



10 



acting with a fluid- wall potential depth of epw = l-OksT [TTl [12] . The bulk was consider 



in z regions away from the walls, at Za = 20a. So Eq. (32) allows the evaluation r^"''^ for 



a layer of width L. The corresponding molecular dynamics r is evaluated by fitting the 
persistence probability curve with Eq. (pi). This fitting was good for all distances z and 
even for different values of L fl] . 

For a typical simple Argon like fluid with a Z}''"''^ of 4.0676 x 10^^ m? / s, for a layer of 
width L = 0.5a = 0.1668 nm, we get a value 0.57002 ps for t^"^^^, while the corresponding 
^buik fj^Qj^ ^YiQ ]\/[j3 simulation is 1.3678 ps. Therefore, in the bulk the MD persistence time 



is more than twice larger that the SE time. From this results, Eqs. (29) and (32) provide a 



procedure to infer about r^^ simply by knowing the local equilibrium particle density inside 



the layer. Rewriting Eq. (32) 



r^bulk bulk t-j* bulk 

SE MD 



12' 



(33) 



1 o^bulk ' 



where the constant D* 

As discussed in [2| , the numerical MD value of r^""^ gives a value of D* = 
1.69512 X 10~^ m'^/s which differs from the expected value of the bulk diffusion constant 
jjbuik ^ This is due to the known difference in the MD and SE average persistence probabili- 
ties and, by definition, of r^^ and t^^^, for the same value of the virtual layer width L. Since 
this is a systematic difference, we conjecture that this relation holds also for layers located 
at Za in the vicinity of an attractive wall 



D{za)r^ 



This can be rearranged as 



XZa] 



T 



bulk 



^ '^MD\^0-)l 



L^ 



12 



r 



bulk 



bulk bulk 
SE 



K 



r. 



iZa) 



r-bulk 



(34) 



''"sBV^aj 



r' 



bulk 



D{Za) 
JJbulk 



(35) 



This relationship is our main result. The left hand side of this relationship corresponds 
to a quantity obtained directly from the persistence probability in the Molecular Dynamics 
simulation, using Eqs. ([s]) and ^. The right hand side of Eq. (35) is a quantity that can 



be evaluated independently from the solution of the Smoluchowski equation, Eq. (29). It 
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FIG. 1: Average persistence time relative to the corresponding bulk value, as a function of layer 



position Za- L = 0.5(7. Dots are the result from the MD simulation 



j-bulk 
MD 



111 . The continuous 



line is the result of scaling ( ^f^^i^ I by the factor 



>ul 
' SE 



D(Za) 
jybulk 



according to Eqs. (13) and (29) 



contains the phenomeno logical diffusion coefficient, introduced in the Langevin equation, 



Eq. (14), in the high friction limit, Eq. (15), to relate the fluctuation and dissipation forces. 



D{z) = |(4fj)^- In Fig. 1 we plot the space average mean exit time Tgj^{za) as a function 
of the position of the virtual layer Za- The continuous line is the result of scaling ( '^^iJu ) 
by the factor j^^ according to Eqs. (13) and (29). The quantity [D{za)Tgj^{za)\ was 
numerically evaluated from Eq 



(|29j), using as input only the equilibrium local particle 
density p{z). As we can see, from Fig. [lithe scaling of the SE and MD predicted by Eq. 



(35) gives a very good numerical agreement. 
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VI. CONCLUSIONS 

Once the relaxation time of the persistence probabihty is identified as the space average of 
the mean first passage time, the evaluation of the quantity [D{za) 'T'sEi^a)] is easily obtained 
in terms of the equilibrium local particle density for a highly anisotropic dense fiuid, next 



to an attractive surface. We have obtained a relationship, Eq. (35), which shows the 
correspondence with the VLMD value, for a fixed width of the absorbing virtual layer. The 
simple scaling factor in terms of the anisotropic diffusion coefficient is physically meaningful, 
since D{z) is introduced in the Langevin equation to account phenomenologically for the 
fiuctuation-dissipation theorem. 
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